## *************************************
##          B3 Loess Regression       
##        Hao Zhang & Ye Zhang
##             2025/7/24
## *************************************

# load packages
library(tidyverse)
library(doBy)

# load data
rm(list = ls())
load("../analysis data/firm analysis.RData")

# check unemployment insurance
city <- summaryBy(insurance_dummy + gdp_pressure + gdpgrowth ~ citycode + year, data = firm, FUN = mean)

pdf("gdp_pressure_loess_unemployment.pdf")
par(mfrow = c(1,1))
par(mar=c(5,5,3,2))
city <- city %>% filter(gdp_pressure.mean < 0.3)
city <- city[order(city$gdp_pressure.mean), ]
lo <- loess(insurance_dummy.mean ~ gdp_pressure.mean + gdpgrowth.mean, city, span = 3)
smoothed10 <- predict(lo)
plot(city$insurance_dummy.mean, x = city$gdp_pressure.mean, main = "Loess Smoothing and Prediction", xlab = "GDP Growth Pressure", ylab = "Unemployment Insurance", cex.lab = 1.2, cex.axis = 1)
lines(smoothed10, x = city$gdp_pressure.mean, col="red")
dev.off()

# check medical insurance
city <- summaryBy(medical_dummy + gdp_pressure + gdpgrowth ~ citycode + year, data = firm, FUN = mean)
city <- city %>% filter(year >= 2004)

pdf("gdp_pressure_loess_medical.pdf")
par(mfrow = c(1,1))
par(mar=c(5,5,3,2))
city <- city %>% filter(gdp_pressure.mean < 0.3)
city <- city[order(city$gdp_pressure.mean), ]
lo <- loess(medical_dummy.mean ~ gdp_pressure.mean + gdpgrowth.mean, city, span = 3)
smoothed10 <- predict(lo)
plot(city$medical_dummy.mean, x = city$gdp_pressure.mean, main = "Loess Smoothing and Prediction", xlab = "GDP Growth Pressure", ylab = "Medical Insurance and Pensions", cex.lab = 1.2, cex.axis = 1)
lines(smoothed10, x = city$gdp_pressure.mean, col="red")
dev.off()
